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Abstract: In regression analysis for deriving scaling laws that occur in various scientific 
disciplines, usually standard regression methods have been applied, of which ordinary least 
squares (OLS) is the most popular. In many situations, the assumptions underlying OLS are 
not fulfilled, and several other approaches have been proposed. However, most techniques 
address only part of the shortcomings of OLS. We here discuss a new and more general 
regression method, which we call geodesic least squares regression (GLS). The method is 
based on minimization of the Rao geodesic distance on a probabilistic manifold. For the case 
of a power law, we demonstrate the robustness of the method on synthetic data in the presence 
of significant uncertainty on both the data and the regression model. We then show good 
performance of the method in an application to a scaling law in magnetic confinement fusion. 

Keywords: regression analysis; information geometry; geodesic distance; scaling laws; 
nuclear fusion 
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1. Introduction 

Regression analysis is an essential instrument for data analysis in numerous branehes of seienee. It is 
used for investigating deterministie relations between variables, for model building and for predietion by 
extrapolation to a previously unseen range of the involved variables. In this paper, we foeus on regression 
analysis applied to the estimation of sealing laws. In various seientifie diseiplines, sueh as astronomy, 
biology, eeology and geology, sealing laws are used to eharaeterize the underlying meehanisms at work 
in the respeetive eomplex systems under study. In general, a sealing law deseribes how a quantity of 
interest y seales when ehanging other quantities Xi, X 2 , ■ ■ ■, xp, on whieh it depends. Sealing laws are 
often expressed in terms of a power law: 

7/ = |3oa;^4^..xP^ (1) 

A erueial property of sueh a power law is seale-invarianee, i.e., when multiplying any of the variables 
Xi by a eonstant a, the power law in Equation (1) essentially remains the same, being multiplied only by 
a eonstant 

In nuelear fusion experiments based on magnetie eonfinement of a hot hydrogen plasma, sealing 
laws are erueial for predieting the performanee of future fusion reaetors, whieh will have a larger size, 
magnetie field, plasma density, etc., eompared to present-day experimental deviees [1]. These sealing 
laws ean be estimated on the basis of datasets from multiple fusion deviees, spanning a signifieant part 
of the parameter spaee. Ordinary least squares regression (OLS) eombined with frequentist theory is the 
statistieal workhorse that is employed for this purpose in the vast majority of oases. However, often, there 
is oonsiderable uneertainty in the experimental data, inoluding the prediotor variables, and in the model 
equations (regression model). However, OLS only deals with uneertainty on the response variables and 
does not oover additional oomplieations, ineluding atypieal observations (outliers), heterosoedastioity, 
oorrelations, non-Gaussian distributions, etc. As sueh, OLS regression is often unsuitable for deriving 
sealing laws [2,3], and many seientifie fields oould benefit greatly from a unified regression methodology 
that is flexible and robust and yet relatively simple to implement. 

In order to be able to handle the eomplieations mentioned above, we have developed a new regression 
method, ealled geodesie least squares regression (GLS). It is based on minimization of the Rao geodesie 
distanee between probability distributions on a manifold equipped with the Lisher metrie. In this paper, 
we briefly introduee the method by means of a simple example involving a power law and Gaussian 
noise. We show the good performanee of the method on synthetie data, introdueing outliers in the first 
ease and studying the effeet of a logarithmie transformation of the data in the seeond ease. Linally, 
we present an applieation to the important sealing eoneerning the power threshold for the transition 
into the high eonfinement regime (H-mode) in nuelear fusion experiments based on magnetie plasma 
eonfinement. The details of the quantities involved in this sealing, their experimental determination and 
the underlying physies are not important for the purpose of this paper. Rather, we here aim at showing 
the performanee of GLS on a ehallenging and heterogeneous real-life dataset. 

The paper is struetured as follows. The method of geodesie least squares regression is deseribed in 
Seetion 2, ineluding a short diseussion on ealeulating geodesie distanees on a Gaussian probabilistie 
manifold, within the framework of information geometry. The next seetion, Seetion 3, briefly introduees 
the database that is used in the subsequent regression experiments, in relation to the sealing law for the 
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H-mode power threshold in fusion plasmas. The experiments involving synthetie data are deseribed in 
Section 4, while the real power threshold scaling is derived in Section 5. Section 6 concludes the paper 
and contains an outlook towards future work related to the methodology. 

2. Geodesic Least Squares Regression 

We start by describing the GLS methodology, which was already introduced in [4,5], but here, 
we go into some more detail. We describe a specific form of GLS regression, and it should be 
stressed that various aspects can be generalized, as will be noted accordingly. Furthermore, several 
elements on which GLS is based are also found in other regression techniques. The strength of GLS 
regression is that it integrates many of these aspects in an elegant way, resulting in a method that is very 
general, flexible and robust. From one point of view, GLS is similar to a class of parameter estimation 
methods that are collectively referred to in the statistics community as minimum distance estimation, 
in that GLS minimizes a distance between a parametric model distribution of the data and an empirical 
distribution [6]. We use the Rao geodesic distance (GD) as a similarity measure, which has the advantage 
that it offers an intuitive geometric interpretation. In addition, there are similarities between GLS and 
the generalized linear model [7]. 

We will consider the case of regression with multiple predictor variables (regressors) and a single 
response variable. For this case, we will show that GLS regression can be regarded as a generalization 
of OLS. However, GLS takes place on a probabilistic manifold, whereas classic OLS operates in a 
flat Euclidean space. Indeed, OLS is based on minimizing the difference, i.e., the Euclidean distance, 
between the predicted and measured values of the response variable. Eikewise, GES is based on 
minimizing the GD between distributions on the probabilistic manifold. Therefore, we start by briefly 
introducing some concepts from information geometry related to distance calculation. 


2.1. Distance in Information Geometry 


In information geometry, a parametric family of probability densities is interpreted as a Riemannian 
differentiable manifold [8]. Each point on the manifold corresponds to a specific probability density 
function (pdf) within the family, and the family parameters represent a coordinate system on the 
manifold. The Eisher information (covariance of the score) provides a unique metric tensor. Eor 
a probability model p({a:m}|{0^}) [9] describing a set {xm} of M variables (m = 1,...,M), 
parameterized by a set {0^} of P parameters {k = 1,..., P), the entries gij of the Eisher information 
matrix are given by (no summation): 


9ij{{Q^}) =E ^lnp({a:™}|{0^})^lnp({a:™}|{0^}) 


del 


= -E 




de^ddl 


lnp({a:m}|{0^}) 


i,j,k = 1, 


P. 


The metric provides the basis for distance measurement between pdfs. Specifically, a geodesic 
curve locally minimizes the distance between two points on the manifold equipped with that metric. 
Through calculus of variations, it can be shown that a geodesic is the solution of the following system of 
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nonlinear second-order ordinary differential equations, known in the language of variational analysis as 
Euler-Lagrange equations [10] and in the present context as geodesic equations: 

p 

+ Y. = 0, r = l,...,P. (2) 


Here, the 0* have been parameterized along the geodesic by t and are the Christoffel symbols of 
the second kind, defined through the metric as: 

■pfc _ \ kr f ^9jr dgir _ dQjj \ 

\dd^ ddi aev ’ 


where denotes the components of the inverse metric. The boundary value problem Equation (2) needs 
to be solved assuming the known values of the coordinates at the boundary points of the geodesic. 

Erom the metric and the solution of the geodesic equations, the length Lg of the geodesic curve 
between two distributions with parameter sets {0]^} and { 02 }, i-s-, the geodesic distance between these 
distributions, may be locally calculated as follows (assuming t runs from zero to one): 



where s represents the arc length. In the framework of information geometry, the geodesic distance 
based on the Eisher metric is often referred to as the Rao geodesic distance (GD). 

Coming back to Equations (2) and (3), it should be noted that closed-form expressions for the GD 
are rarely available. On the other hand, provided the Eisher metric can be calculated relatively easily, 
the framework of information geometry is very useful, since straightforward approximations of the 
geodesic curves can be found in a geometrically intuitive way [11]. This intuitive approach by means of 
geometry is an important and attractive aspect of the theory, as it provides enhanced insight into various 
concepts and algorithms in probability theory and statistics [12]. This is also the case for GLS, as will be 
demonstrated below. Eurthermore, as far as the GD is concerned, visualization of geodesics may guide 
controlled approximations to the geodesic paths and geodesic distances [11]. 

Besides the attractive feature of providing intuitive geometrical insight into problems involving 
similarity measurement between probability distributions, the GD has several more advantages over 
other similarity measures for distributions. Eirst, it is a distance measure (a metric) in the strict sense of 
the word. As a result, it is symmetric in its arguments, a desirable property for measuring the similarity 
between two given states of information in terms of probability distributions. In addition, it obeys 
the triangle inequality, yielding various practical advantages, for instance in the field of data retrieval 
from large databases [13]. Eurthermore, closed-form expressions may be available for the GD, or its 
approximation, for various families of distributions where no such analytic form has been found in 
the case of, for instance, the Kullback-Eeibler divergence (KED) [11]. Einally, there is considerable 
experimental evidence suggesting that the GD in general is a more effective similarity measure between 
distributions than the KED (see [11] and the references therein). We note that for distributions that lie 
infinitesimally close on the probabilistic manifold, it can be proven that the Kullback-Eeibler divergence 
equals half of the squared geodesic distance between the distributions (see, e.g., [14]). Hence, in such a 
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case, the KLD and GD yield similar results, but in general, they are quite different measures of similarity 
between distributions. 


2.2. Geodesics for the Univariate Normal Distribution 


In this paper, we diseuss applieations that are based on a univariate normal distribution c^), 
parameterized by its mean p. and standard deviation cr. In this ease, an analytie expression for the 
Fisher-Rao metrie is available. It turns out to be the familiar Poineare metrie, whieh, when applied to a 
half-plane, is a well-known model for hyperbolie geometry that has eonstant negative sealar eurvature. 
The line-element is given by [15,16]: 




dp.^ da^ 
— -f 2- 

rr2 ^ rr2 


(4) 


As an illustration, the Poineare half-plane is pietured in Figure la, together with two geodesies 
between, on the one hand, the points pi = M (4,1.2^) and p 2 = M (16,1.5^) and, on the other hand, 
P 3 = AT (4,4.0^) and p 4 = Af (16, 5.0^). The eorresponding normal density funetions are drawn in 
Figure lb, as well as a number of densities assoeiated with some intermediate points on eaeh geodesie. 
As a further illustration. Figure le shows one blade of a partieular pseudosphere, namely the traetroid, 
whieh is loeally isometrie to the Poineare half-plane and the univariate normal manifold for ff > 1, 
with periodieity in p. In order to better visualize the geodesies, a resealed version of the traetroid 
is shown in Figure Id. This surfaee has a longer period in the p-direetion. However, it should be 
kept in mind that only the visualization in Figure le ean be used to measure absolute distanees on the 
surfaee, for in Figure Id, the pietured geodesies are no longer the shortest eurves between the points 
in question. It is elear that the geodesies on the Gaussian manifold are different from straight lines 
in the Euelidean spaee, wherein the manifold has been immersed. The shape of the geodesies ean be 
made intuitively elear by noting that they always pass through a region of inereased standard deviation 
relative to that of the boundary points. This provides the shortest route, as ean be seen from the line 
element Equation (4). Interestingly, similar arguments will be shown to enable a deeper insight into the 
operation of GLS regression. We further note that various alternative models exist to visualize hyperbolie 
geometry; see, e.g., [17]. 

A elosed-form expression is available for the GD on the normal manifold, permitting fast evaluation. 
Indeed, for two univariate normal distributions pi (x| p-i, erf) and p 2 (x| p. 2 , erf), the GD is given by [16]: 


1 + 6 


GD(pi,P 2) = V^ln-- = 2A/2tanh ^6, 6 = 


1 - 6 


(p,i - p2)^ + 2(0-1 - 0'2)' 
_(pi - P2)2 + 2(0-1 + o-2)2_ 


1/2 


(5) 


Furthermore, sinee the injeetivity radius of the hyperbolie plane is infinite, the geodesies are globally 
length-minimizing [10]. 
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Figure 1. (a) Illustration of the Poincare half-plane with several half-circle geodesics in the 
background, together with the geodesic between the points pi and p 2 and between ps and p 4 , 
defined in the main text, (b) Probability densities corresponding to the points pi, p 2 , pa and 
P 4 indicated in (a). The densities associated with some intermediate points on the geodesics 
between pi and p 2 and between pa and p 4 are also drawn, (c) Rendering of one blade of the 
tractroid, again with the two geodesics superimposed. The parallels of the tractroid are lines 
of constant standard deviation cr, while the meridians (the tractrices) are lines of constant 
mean p.. This representation of the normal manifold is periodic in the p.-direction, and a 
rescaled version (longer period along p.) is shown in (d). 


2.3. Geodesic Least Squares Methodology 

GLS starts from the premise that the probability distribution underlying experimental measurements 
is the fundamental object resulting from the measurement. As such, GLS does not perform regression 
based on data points in a Euclidean space, but rather operates on probability distributions lying on 
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a probabilistic manifold. This introduces additional flexibility that renders the method robust in the 
presenee of large uneertainties, as will be demonstrated in the experiments. 

Briefly, the idea is to eonsider two different proposals for the distribution of the response (dependent) 
variable y, eonditional on the predietor variables. On the one hand, there is the distribution that one would 
expeet if all assumptions were eorreet regarding the deterministie eomponent of the regression model 
(regression funetion) and the stoehastie eomponent. We eall this the modeled distribution. On the other 
hand, we try to eapture the eonditional distribution of y by relying less on the model assumptions, but 
direetly on the measurements of y. For this, we will use the term observed distribution. It is in this sense 
that GLS is similar to minimum distanee estimation (MDE), where the Hellinger distanee is a popular 
similarity measure [18]. This was first applied to regression in [19], but there are several differenees with 
GLS. First and foremost, GLS ealeulates the geodesie distanee between eaeh individual pair of modeled 
and observed distributions of the response variable, eorresponding to an individual measurement point. 
As sueh, eaeh individual data point aequires the status of a probability distribution in its own right. 
Consequently, GLS performs regression between probability distributions on a probabilistie manifold. 
In eontrast, MDE usually eonsiders a distanee between a kernel density estimate of the distribution of 
residuals, on the one hand, and the parametrie model, on the other hand, based on the entire data sample. 
Seeondly, we explieitly model all parameters of the modeled distribution, whieh is similar to the ideas 
behind the link funetion in the generalized linear model (GLM) [7]. In the present work, this will be 
aeeomplished by explieitly modeling both the mean and standard deviation of the Gaussian modeled 
distribution. Additionally, a final differenee is that we use the Rao geodesie distanee as a similarity 
measure. 

As a simple example that we will use also in the experiments, eonsider a linear relation p = (3^ 
between a single predietor variable £, and a response variable p, with |3 a eonstant. In aeeordanee 
with the diseussion above, we explieitly wish to allow for the ehallenging ease of uneertainty on the 
predietor variable £,. Therefore, we assume that, in reality, N samples of a stoehastie (noisy) variable x 
are observed, together with N samples of a stoehastie response variable y. We take the simple ease of 
normally distributed (Gaussian) noise: 

y = ^ + ey = ^E, + ey, (O, aj) , 

( 6 ) 

x = E, + e^, ~ AC (O, 0-^) . 

The observations Xn (n = 1,..., N) are taken as mutually independent and so are the ?/„. and 

Oy are assumed to be known, and in this example, they are taken eonstant for all measurements, i.e., 
we have homoseedastieity. However, we will also eonsider heteroseedastieity later on. Aeeording to the 
regression model, eonditionally on x„, eaeh measurement is drawn from a normal distribution: 

Pmod{y\xn) = AC (|3x„, ) where (7) 

with the subseript “mod” referring to the modeled distribution. In our simple example. Equation (7) 
follows from standard Gaussian error propagation rules. However, for nonlinear regression laws, 
the eonditional distribution for y has to be obtained by marginalizing the unknown true values 
Nevertheless, the Gaussian error propagation laws may be used in the nonlinear ease as well, 
to approximate the eonditional distribution p{y\xn) by a normal distribution, as will be shown 
in the experiments. 
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We next ehoose a speeifie form of the observed distribution eorresponding to eaeh realization of the 
variable y, eonditional on the observations, i.e., PohsiylVn)- In this example, we take again the normal 
distribution, but eentered on eaeh data point: Af (j/„, cr^bs)’ where cTobs is to be estimated from the data. In 
the eontext of the GLM, this is known as the saturated model. The extra parameter O'obs gives the method 
added flexibility, sinee it is not a priori required to equal cTmod- As a result, GLS is less sensitive to 
ineorreet model assumptions. Note that in this example, we have ehosen the observed distribution from 
the same model (Gaussian) as the modeled distribution. Furthermore, O'mod is taken as a fixed value for 
all measurements and so is O'obs- These assumptions ean of eourse be relaxed, leading to a more general 
method. However, the transition from OLS to GLS is best explained by means of a Gaussian observed 
distribution, whieh, in addition, offers eomputational advantages, sinee the expression for the GD has a 
elosed form; see Equation (5). 

GLS now proeeeds by minimizing the total GD between, on the one hand, the joint observed 
distribution of the N realizations of the variable y and, on the other hand, the joint modeled distribution. 
Thanks to the independenee assumption in this example, we ean write this in terms of produets of the 
eorresponding marginal distributions: 




argmin GD 

(3 GM 


■ N N 

J_ J_ Pohs{y\yn) 1 PmoA{,y\Xn) 
_n=l n=l 


N 

argmin ^ GD^ [Pobs(2/|2/n),Pmod(l/kn)] • 
p6R 


( 8 ) 


The last equality entails a eonsiderable simplifieation, owing to the property that the squared GD 
between produets of distributions ean be written as the sum of squared GDs between the eorresponding 
faetors [16]. Henee, the optimization proeedure involves matehing not only yn with hxn, but also 
cr^bs with Note that the parameter |3 oeeurs both in the mean and the varianee of the 

modeled distribution. Ineidentally, foreing a^bs = cr^ + |3^cr^ would take us baek to standard maximum 
likelihood estimation, for the Rao GD between the two Gaussians pobs and Pmod with means yn and 
bxn, respeetively, but with identieal standard deviations (fixed along the geodesie path), is preeisely the 
Mahalanobis distanee [20] : 


GD(pobs) Pmod) 


\yn - hXn\ 


if o'is 


= S + P 


V!. 


We note that the GLS seheme addresses many of the diffieulties with elassie OLS regression. Lirst, 
GLS explieitly allows uneertainty on the predietor variables, and it is not restrieted to normal or 
symmetrie noise distributions, nor does it neeessarily assume homoseedastieity. In addition, eorrelations 
among variables and among observations ean be built into the stoehastie eomponent of the regression 
model. Lurthermore, GLS ean operate with any (nonlinear) regression funetion. Moreover, it will be 
shown in the experiments that GLS is relatively insensitive to uneertainties in both the stoehastie and 
deterministie eomponents of the regression model. The same quality renders the method also robust 
against outliers. 

In the experiments below, we employed a elassie aetive-set algorithm to earry out the 
optimization [21]. Lurthermore, presently, the GLS method does not direetly offer eonfidenee (or 
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credible) intervals on the estimated quantities. Future work will address this issue in more detail, but 
for now, error estimates were derived by Monte Carlo sampling in the ease of the numerieal simulations 
(Seetion 4) and by bootstrapping in the ease of the real data (Seetion 5) [22]. The bootstrapping involved 
ereating, from the measured dataset, a large number of artifieial datasets of the same size, by resampling 
with replaeement. The regression analysis was then earned out on eaeh of the datasets, and the mean 
and standard deviation, over all datasets, of eaeh estimated regression parameter and of the predieted 
quantities were used as estimates of the parameter or predietion value and its error bar, respeetively. 
This seheme typieally results in rather eonservative error bars, whieh eould possibly be narrowed down 
using more sophistieated methods. 

3. The L-H Power Threshold and Database 

We now provide some baekground information regarding the main regression applieation that will be 
treated in the experiments with synthetie and real data. It eoneems one of the most important sealing 
relations in fusion seienee based on magnetie plasma eonfinement, related to the threshold Pthr for the 
heating power that is required for the plasma to make the transition into a desired regime of high energy 
eonfinement (H-mode) in the next-step fusion deviee ITER (International Thermonuelear Experimental 
Reaetor) [1,23,24]. To a good approximation, this so-ealled L-H (or H-mode) power threshold depends 
on the eleetron density in the plasma he (in 10^° m“^), the main magnetie field Bt (in tesla (T)) and the 
total surfaee area S of the eonfined plasma (in m^). This is usually expressed by means of the following 
sealing relation: 

Fthr = (9) 

To estimate the eoeffieients in this relation, we employed data from eight fusion deviees of the 
tokamak type (ASDEX, AUG, CMOD, DIII-D, JET, JET-2M, JT-60U, PBXM) in the International 
Tokamak Physies Aetivity (ITPA) multi-maehine database for the L-H power threshold (subset 
IAEA02 [23,25-27]). This yields 616 measurement sets of power, density, magnetie field and surfaee 
area, eaeh set obtained during a brief time of plasma operation under stationary eonditions in one of the 
eight deviees involved in the study [28]. 

The ITPA database eontains some information regarding the error bars on the measurements, 
speeifieally relative errors expressed as pereentages. This is important for our purposes, beeause we 
need the error bars to ealeulate ffmod- Unfortunately, the error estimates are not available in some eases, 
and if they are, the preeise definition of the error bars is not always elear. Usually, an error bar in the 
database represents an estimate by the experimentalist of the typieal range in whieh the “true” quantity 
ean be expeeted to lie, where the uneertainty is assumed to be eaused by both stoehastie and systematie 
effeets. Moreover, it is diffieult to assess the probability that is eovered by the stoehastie eomponent of 
the errors mentioned in the database. Sinee a detailed investigation of the uneertainty of the threshold 
data is beyond the seope of the present paper, we will assume that the error bars pertain to a stoehastie 
uneertainty eorresponding to a single standard deviation of a Gaussian distribution. Eor some derived 
quantities, the error bars had to be ealeulated from the uneertainty on more fundamental measurements. 
In those eases, we employed Gaussian error propagation rules to estimate the standard deviation on the 
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derived quantities. For the ease of the global H-mode eonfinement database, this strategy has been shown 
to provide reasonable information on the aetual error bars [29] . 

It is important to mention that the main souree of uneertainty in the data used for power threshold 
sealing, when eompared to the predietions of a simple power law regression model, is not expeeted to be 
due to the measurement uneertainty on the individual variables. There are far larger sourees owing to the 
variability of the experiments that produeed the data. To estimate the variability of eaeh of the physieal 
quantities with respect to the scaling law, we performed the following calculation. First, the nonlinear 
scaling law was estimated using OLS, as explained in Section 5.2. Then, for a specific variable z (one 
of the predictor variables or the dependent variable) and for each data point, the relative difference was 
computed between the 2 ;-value of the data point itself and the 2 ;-value of the projection of the data point on 
the hypersurface given by the scaling law, keeping the values of the other variables fixed. This difference 
can be interpreted as the deviation of the point from the theoretical scaling law, assuming the deviation is 
solely due to the variability of the variable z. Finally, the standard deviation of these relative differences 
was taken, and the procedure was repeated for every predictor variable and the dependent variable. The 
resulting standard deviations can be interpreted as upper bounds of the relative variability of each of the 
variables around their ‘theoretical’ values given by the scaling law. This way, for Ug, Bt, S and Pthr, we 
obtained 39%, 31%, 28% and 38%, respectively. These levels are clearly much larger than the relative 
uncertainties due to measurement error alone. Indeed, the typical measurement error bars quoted in the 
ITPA database, on average, over all devices, are estimated at 4% for h-e, 1% for Bt, 3% for S and 15% 
for Pthr [25,26]. 


4. Numerical Simulations 

We next demonstrate some of the potential of the GLS regression scheme by means of a number 
of experiments with synthetically-generated data. We treat two particular cases of deviation from the 
model according to which the data were created and show that, in comparison with a number of standard 
regression techniques, GLS yields the most accurate results across all experiments. The first case 
concerns the effect of outliers, while in the second case, the influence of a logarithmic transformation is 
studied. In each case, we start with a very simple experiment that is easily reproduced, using a single 
predictor variable, providing some intuitive feeling regarding the performance of the method. We then 
proceed to a more elaborate test, still based on partly synthetic data, but using a regression challenge 
similar to that used in the real-world experiment for scaling of the L-H power threshold in fusion plasmas, 
presented in Section 5. 

4.1. Effect of Outliers 

The robustness of minimum distance estimators to outliers in the data was noted in the classic 
literature of minimum distance estimation [18]. We now show that this is a quality also enjoyed by 
GLS regression. 
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4.1.1. Single Predictor Variable 

We first concentrate on estimating the slope of a regression line with a single predictor variable. 
To this end, a dataset was generated consisting of ten points labeled by coordinates and 
r|„ (n = 1,..., 10), with the £,„ chosen unevenly between zero and 50 and r\n = (3£,n, taking |3 = 3. 
Then, Gaussian noise was added to all coordinates according to Equation (6), with cSy = 2.0 and 
= 0.5, resulting in values for the predictor variable and for the response variable. Finally, 
one outlier was created by multiplying the value of yk by a factor distributed uniformly between 1.5 and 
2.5, with k chosen uniformly among the indices 8, 9 and 10. 

We next estimated (3 by means of GLS and compared the estimates with those obtained by OLS, 
maximum a posteriori (MAP) using the model in Equation (7) for the likelihood and an uninformative 
prior [30], total least squares (TES), which is a typical errors-in-variables technique [31], and a robust 
method (ROB) based on iteratively reweighted least squares (bisquare weighting) [32], included in the 
MATEAB Statistics toolbox [33]. It should be noted that MAP takes into account the error bars on the 
predictor variables. In all cases, we assumed knowledge of the values of and Gy. In order to get an 
idea of the variability of the estimates, Monte Carlo sampling of the data-generating distributions was 
performed, and the estimation was carried out 100 times. 

The results are given in Table 1, mentioning the sample average and standard deviation of the 
estimates |3 over the 100 runs for each of the methods. GES is seen to perform very well and similar to 
the robust method ROB, but the other techniques yield considerably worse results. The average estimate 
of CTobs was 5.43 with a standard deviation of 0.24. On the other hand, the modeled value of the standard 
deviation in the conditional distribution for y was ffmod = y/o|”+”9o^ = 2.5. Hence, GES succeeds 
in ignoring the outlier by increasing the estimated variability of the data. Put differently, the effect of 
the outlier is, in a sense, to increase the overall variability of the data, which GES takes into account 
by increasing the observed standard deviation of the data (aobs) with respect to the standard deviation 
predicted by the model (cTmod)- 

Table 1. Monte Carlo estimates of the mean and standard deviation for the slope parameter 
in linear regression with errors on both variables and one outlier. GES, geodesic least squares 
regression; TES, total least squares; ROB, robust method. 


Original 

GLS 

OLS 

MAP 

TLS 

ROB 

o 

o 

CO 

II 

3.031 ± 0.035 

3.68 ± 0.29 

3.83 ± 0.36 

4.6 ± 1.0 

2.992 ± 0.041 


As mentioned before, this result can also be understood in terms of the pseudosphere as a geometrical 
model for the normal distribution. To see this, we refer to Figure 2, where several sets of points 
(distributions) are drawn on a portion of the surface of the pseudosphere for one particular dataset 
generated as described above. First, the modeled distributions are plotted with their means |3x„ 
(see Equation (7)) and standard deviations Cmod = 2.5, using the average estimate p = 3.031 obtained 
by GES. These are the green points on the surface, and they lie on a parallel, since they all correspond 
to Gaussians with the same standard deviation ffmod- In this particular dataset, the index of the outlier 
was k = 10, so the point |3xio is indicated individually. Obviously, according to the model, no outlier is 
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expected, so the modeled distribution corresponding to k = 10, which is the green point just mentioned, 
lies close to the other predicted points (distributions). Next, we plot the observed distributions with their 
means Pn and standard deviations ffobs (for this dataset estimated at dobs = 5.43). These are the blue 
points, lying at a constant standard deviation Oobs. which is higher than ffmod (5.43 > 2.5). The outlier 
2/10 can clearly be observed, and being an outlier, it lies relatively far away from the rest of the blue 
points (observed distributions). Now suppose that, like MAP, GLS would not be able to increase ffobs 
relative to cTmod in order to accommodate the outlier. Then, the observed distributions would have the 
same observed means (the measured values pn), but they would have the standard deviation predicted 
by the model. Hence, they would lie on the parallel corresponding to O'mod, just like the green points. 
We have plotted these fictitious distributions as the red points at the level of cTmod, and they are labeled 
p. Again, the outlier (labeled ^lo) can be seen, but it seems to lie further away from the other red points 
(the points pn) compared to the actually observed situation, i.e., the distance from piQ to the other p^ 
(blue points). At least this is the case when using (visually) the Euclidean distance in the embedding 
Euclidean space. We can verify that this is indeed so by using the proper geodesic distance on the 
surface: overall, the blue points lie closer together (including the outlier) than the red points. Now, 
in fact, GES aims at minimizing the distance between each green point (modeled distribution) and its 
corresponding blue point (observed distribution), so as far as the outlier is concerned, we should really be 
looking at the geodesic between the point (pxio, o'mod) and the point (//lo, cTobs)- The geodesic (labeled 
“Geoi”) between these points is also drawn on the surface, and again, we compare this to the fictitious 
situation, represented by the geodesic (labeled “Geo 2 ”) between (^xio, cTmod) and (^lo, o'mod)- Indeed, 
again, we see that the geodesic Geoi is shorter than Geo2. Therefore, by increasing cTobs relative to amod, 
the outlier is not so much an outlier anymore, as measured on the pseudosphere!When calculating the 
GD, one finds 2.4 for Geoi and 2.8 for Geo2. Therefore, GES obtains a lower value of the objective 
function (sum of squared geodesic distances) if it increases ffobs with respect to ffmod- Of course, there 
is a limit to this: GES cannot continue raising ffobs indefinitely, trying to mitigate the distorting effect 
of the outlier, for then, the other points would get a too high observed standard deviation, which is not 
supported by the data. The image that we see in Eigure 2 is the best compromise that GES could find. In 
fact, we note that, in the case we suspect that pi^ could be an outlier, it may very well be worthwhile to 
introduce two parameters to describe the observed standard deviation: one for the nine points that seem 
to follow the model and one to take care of the outlier. This would be a very straightforward extension of 
the method, and we explore this to some extent when using data from the ITPA database below. There, 
we assign a separate parameter to describe the observed standard deviation of all data coming from a 
specific tokamak, hence defining an individual parameter for each machine. 


Entropy 2015, 17 


13 



Figure 2. A portion of the pseudosphere together with the regression results on synthetie 
data with an outlier, as described in the main text. 


4.1.2. Multiple Predictor Variables 


In this experiment, a regression problem with multiple predictor variables and a power law is studied. 
The deterministic part of the regression model is based on the real-world problem for the L-H power 
threshold in fusion plasmas, which we are going to consider in Section 5. Furthermore, the values of 
the predictor variables are taken from the same international power threshold database, and values of the 
response variable are synthetically generated from this. 

More specifically, the dataset for this experiment was created as follows. First, an artificial linear 
regression law was put forward for a variable p, depending on the predictor variables he, and S, 
which were introduced in the context of the power threshold scaling law in Section 3 [34]. In particular, 
we generated a number of realizations of the variable p from the following prescription: 

^ = 1^0 + Pihe + ^2Bt + Pas'. (10) 


This was considered as the “true” relation between the predictor and response variables, where, as 
mentioned above, the values of the predictor variables were chosen to be exactly those from the ITPA 
database, which are normally used in the real power threshold scaling law. A whole range of datasets 
was created using the following values of the coefficients po. Pi> P 2 and P 3 : 


Po = 1,1.1,..., 20, 

Pi, P 2 , Ps = 0.1, 0.2,..., 2. 


( 11 ) 


Thus, for each combination of values of Po, Pi, P 2 and ps, all 616 values of r\ were calculated 
according to Equation (10), based on the values of he, Bt and S from the ITPA database. The range 
of coefficient values in Equation (11) was chosen to be representative for the values that are typically 
obtained from a regression analysis on the true scaling law (see Section 5). The exception is Po, for 
which the range was chosen of roughly the same order as p — po (much smaller values of Po would not 
be estimable in comparison with p — Po). 
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Next, Gaussian noise was added to both the predietor and response variables. The noise level was 
ehosen according to the typical relative measurement errors in the ITPA database, i.e., 4% for ffe, 
resulting in a variable Xi, 1% for (variable X 2 ), 3% for S (variable Xs) and 15% for the dependent 
variable (variable y, which is Pthr in the real-world regression problem). It should be stressed that, in the 
light of our comments in Section 3 regarding the variability of the predictor variables, these are rather 
low noise levels. We further note that fixed relative noise levels lead to a different standard deviation for 
each measurement (heteroscedasticity). 

Furthermore, in this experiment studying the effect of atypical observations, 10 outliers were created 
in each of the datasets. In particular, from the total of 616 points in each dataset, 10 points were randomly 
chosen, and the associated value of y was multiplied with a factor F, where F was distributed uniformly 
between 1.5 and 2.5. For each combination of coefficient values (3* (i = 0,..., 3) taken from Equation 
(11), 10 datasets were realized, each time performing the sampling of noise and outliers. 

Finally, the regression analysis was carried out for every dataset, and for each choice of the |3j, the 
obtained estimates |3j were defined as the average over the 10 data realizations. Next, histograms were 
created based on these averages for the estimated coefficients, specifically the normalized histograms of 
the relative difference (|3j — |3j)/(3i (i = 0,..., 3), expressed as a percentage, between the true value (3* 
and the estimated value |3j of each regression parameter. The histograms of these percentage errors are 
shown in Figure 3. In order to avoid a cluttered figure, the results of OLS, MAP and GLS are plotted in 
one panel and those of TLS and ROB in another. 
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Figure 3. (a) Histograms of the relative error in estimating the regression coefficients |3i by 
means of OLS, MAP and GLS for a linear regression problem with outliers. Horizontal axes 
represent the error in percent and vertical axes probability, normalized to one. (b) Similar, 
for TLS and ROB. 
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From these histograms, it is elear that, for eaeh parameter, GLS performs mueh better than OLS, 
MAP and TLS, with the latter failing eompletely. In ease of GLS, the vast majority of relative errors 
is of the order of a few pereent and eertainly smaller than 20%. Overall, the most diffieult to estimate 
parameter turns out to be (3i, whieh is assoeiated with he. The robust estimation teehnique in MATLAB 
also delivers good results (in faet, not mueh worse than GLS), as it is designed to eope with outliers. 
However, we will see that in the next experiment ROB does not perform well at all. 

4.2. Effect of Logarithmic Transformation 

We next tested the effeet of a logarithmie transformation, whieh is often used to transform a power-law 
regression model into a linear form. However, the logarithm alters the data distribution, whieh may lead 
to misguided inferenees from OLS [2,3]. Therefore, the flexibility offered by GLS is expeeted to be 
benefieial in this ease, as it allows the observed distribution to deviate from the modeled distribution. 

4.2.1. Single Predietor Variable 

Again, we first performed a simple regression experiment involving a single predietor variable, with 
a power law deterministie model and additive Gaussian noise on all variables. In aeeordanee with 
the typieal situation of fitting fusion sealing laws to multi-maehine data, the noise standard deviation 
was taken proportional to the simulated measurements, eorresponding to a given set of relative error 
bars. As a result, in the logarithmie spaee the distributions were only approximately Gaussian, with the 
standard deviation given by the eonstant relative error on the original measurement (homoseedastieity). 
Ten points were ehosen with predietor values unevenly spread between zero and 60. A power law 
was proposed to relate the unobserved and r\n' 

r\n = ^o^t\ n = l,...,10. 

Then, Gaussian noise was added to the £,„ and ri„, eorresponding to a substantial relative error of 
40%. We finally took the natural logarithm of all observed values and ?/„, enabling applieation of the 
same linear regression methods that were used in the previous experiment. In this partieular experiment, 
we ehose (3o = 0.8 and (3i = 1.4, but we found that other values yield similar eonelusions. Again, 
100 data replieations were generated, allowing ealeulation of Monte Carlo averages. 

The averages and standard deviations over all 100 runs are given in Table 2. Again, the results show 
that GLS is robust against the flawed model assumptions, now performing similar to TLS. 

Table 2. Monte Carlo estimates of the mean and standard deviation for the parameters in a 
log-linear regression experiment with proportional additive noise on both variables. 


Parameter 

Original 

GLS 

OLS 

MAP 

TLS 

ROB 

Po 

0.80 

0.94 ± 0.47 

2.2 ± 2.3 

3.0 ± 1.7 

0.99 ± 0.70 

2.72 ± 0.77 

Pi 

1.40 

1.39 ± 0.11 

1.19 ± 0.16 

1.08 ± 0.26 

1.41 ± 0.14 

1.17 ± 0.11 
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4.2.2. Multiple Predictor Variables 


In the last experiment with synthetic data, we studied the effect of a logarithmic transformation in a 
similar problem as the one described in Section 4.1.2, but in the case of a power law. In particular, the 
variable p was calculated for the same range of values of the parameters (3* as given in Equation (11), 
but now according to a power law: 

p = 

Then, Gaussian noise was added to all variables. However, when applying the relatively low noise 
levels used in Section 4.1.2, no significant differences were observed in the performance of GLS and 
MAP. Therefore, the noise levels for the predictor variables were augmented to 20% for Pe (variable xi), 
5% for Bt (variable X 2 ) and 15% for S (variable X 3 ). The level for Pthr was kept at 15%, as before. 
This is still well within the maximum variability range that can be expected for the predictor variables in 
the ITPA database, as discussed in Section 3. 



Figure 4. Histograms of the relative error in estimating the regression coefficients |3j by 
means of OLS, MAP and GLS for a power-law regression problem after a logarithmic 
transformation. Horizontal axes represent the error in percent and vertical axes probability, 
normalized to one. (b) Similar, for TLS and ROB. 


After adding the noise, all data were transformed to the logarithmic domain, and 10 datasets were 
generated for each combination of regression coefficients. Subsequently, linear regression analysis was 
applied to each of the log-transformed datasets. The coefficient estimates, defined as the average over the 
10 replications, were then compared among the various regression methods, as shown in Figure 4. Again, 
the normalized histograms of the relative error on the estimated parameters are displayed, showing 
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the consistently better performance of GLS over all other methods tested, including TLS and ROB. 
For GLS, the errors on (3o and |3i are the largest, compared to those on (32 and (Sa, but the majority is 
still below 20%. As for (3o, the slightly inferior performance of GLS relative to the results with outliers in 
Section 4.1.2 is simply due to the fact that log (Sq for the lowest values of (3o is negligibly small compared 
to logri - log (3o. 

5. Power Threshold Scaling 

We finally come to the application of power threshold scaling using real-world data from the ITPA 
database for all variables, including the response variable Pthr- We start with log-linear regression and 
then apply nonlinear regression analysis. Next, we perform a simple analysis of the influence of the error 
bars on the estimation results, and we finally provide a discussion of the results in this section. 

5.1. Linear Scaling 

We first followed the standard practice of transforming to the logarithmic scale to estimate the 
coefficients (3o, (3i, (32 and (33 in Equation (9) via linear regression. In the GLS method, we introduced 
additional parameters Cobs,a N^), one for each of the Nt = 8 tokamaks contributing data to 

the scaling. That is, if a certain data point with index n originated from tokamak a, then in term n of 
the objective function in Equation (8), an observed distribution was used, parameterized by means of the 
O'obs,a corresponding to that machine. The Cobs,a serve a similar purpose as the parameter ffobs defined 
above, except that they describe the observed standard deviations of the logarithmic power threshold. 
This, of course, corresponds to the relative errors on the power threshold itself. To calculate ffmod for 
each data point, we used the relative measurement error bars quoted in the database (typically 4% for n^, 
1% for Bt, 3% for S and 15% for Pthr)- Considering the discussion in Section 3 regarding other sources 
of uncertainty, it is clear that the Cobs.o will need to take into account other, “unexpected” uncertainty 
sources, hence increasing the flexibility of the method. 

In this analysis, we compared GLS only with OLS and the powerful MAP method. The results on the 
IAEA02 data are given in Table 3. The predictions for ITER are also shown, for two typical densities 
(0.5 and 1.0 X 10^° m ^). All estimates are accompanied by their 95% credible intervals obtained from 
100 bootstrap samples (artificial datasets). We stress that this notion of a credible interval corresponds to 
the standard Bayesian definition of an interval wherein the true value of a stochastic variable is assumed 
to lie with a certain probability (e.g., 0.95). 

The estimates by GES of the parameters Cobs.a (observed standard deviation on log Pthr), for each of 
the devices contributing to the IAEA02 data, were expressed as a relative error on the bootstrap-averaged 
Pthr- These relative errors and their credible intervals are given in Table 4. The relative error on the power 
threshold lies around 15% to 30% for the various machines, except for ASDEX, where the uncertainty 
reaches a higher level of about 40%. On average, this yields an estimated error of 24.2% for Pthr, which 
is quite somewhat higher than the average of 15% mentioned in the database, although still considerably 
lower than the upper bound of 38%, as calculated in Section 3. Again, this is an indication of additional 
sources of uncertainty, on top of mere measurement error, causing the data points to deviate from the 
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proposed regression model, as diseussed already in Seetion 3. That extra uneertainty is eaptured by the 
GLS method. 


Table 3. Estimates of regression parameters and predietions for ITER in log-transformed 
linear sealing of the H-mode threshold power using the IAEA02 dataset. The bootstrap 
averages are given, as well as the 95% eredible intervals (Cl). 










Method 


Po 

Pi 

P2 

Pa 

Pthr,0.5 (MW) 

Pthr, 1.0 (MW) 

OLS 

Average 

0.0507 

0.485 

0.873 

0.843 

38.0 

53.2 


Cl 

±0.0060 

±0.073 

±0.061 

±0.041 

±4.4 

±8.0 

MAP 

Average 

0.0449 

0.567 

0.867 

0.901 

45.6 

67.6 


Cl 

±0.0051 

±0.078 

±0.069 

±0.039 

±5.0 

±9.6 

GLS 

Average 

0.0426 

0.660 

0.795 

0.946 

48.3 

76.4 


Cl 

±0.0042 

±0.069 

±0.059 

±0.034 

±4.7 

±9.8 


Table 4. Estimates of the observed standard deviations O'obs.a of tho logarithmie power 
threshold, expressed as pereentage errors on Pthr itself, for the tokamaks eontributing to 
the IAEA02 dataset, obtained using log-transformed linear sealing. The bootstrap averages 
are given, as well as the 95% eredible intervals (Cl). 



ASDEX 

AUG 

CMOD 

DIII-D 

JET 

JFT-2M 

JT-60U 

PBXM 

Average (%) 

41.8 

23.0 

22.0 

15.7 

24.6 

15.9 

22.8 

27.6 

Cl (%) 

±5.3 

±1.4 

±1.1 

±1.8 

±2.0 

±1.2 

±2.3 

±2.9 


5.2. Nonlinear Scaling 


Next, we show the results of nonlinear regression in the original data spaee, i.e., without logarithmie 
transformation. Whereas this prevents an analytie solution using OES, the advantage is that the 
distribution of the data is left undistorted [2,3], while the implementation of OES, MAP and GES is 
not signifieantly more eomplex. Indeed, the distribution of the right-hand side in Equation (9) ean be 
approximated by a Gaussian with mean p^od = Po and standard deviation Cmod, given by: 


2 2 2 
^mod ^^thr l^mod 




Henee, the modeled standard deviations depend on the measurements (heteroseedastieity). 
Nevertheless, in defining the observed standard deviations crobs,a> we introdueed an approximation 
assuming eonstant error bars for all measurements from a single maehine. This assumption may be 
relaxed in the future. 

The results of the sealings and predietions are presented in Tables 5 and 6. We eompared GES with 
OES and MAP using uniform priors. It may be possible to derive even less informative priors for MAP, as 
was done in the log-linear ease in Seetion 5.1 (and see [30,35]), but this was not pursued here. Moreover, 
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even in the log-linear analysis, we observed only a marginal differenee between the results under various 
ehoiees of priors. 


Table 5. Estimates of regression parameters and predietions for ITER in power-law sealing 
on the original seale of the H-mode threshold power using the IAEA02 dataset. The bootstrap 
averages are given, as well as the 95% eredible intervals (Cl). 










Method 


Po 

Pi 

P2 

P3 

Pthr.0.5 (MW) 

Pthr,1.0 (MW) 

OLS 

Average 

0.0274 

0.773 

0.96 

1.038 

69 

118 


Cl 

±0.0083 

±0.090 

±0.10 

±0.071 

±15 

±32 

MAP 

Average 

0.0425 

0.643 

0.788 

0.933 

44.2 

69.1 


Cl 

±0.0041 

±0.074 

±0.079 

±0.034 

±3.8 

±8.2 

GLS 

Average 

0.0397 

0.715 

0.751 

0.984 

51.6 

84.7 


Cl 

±0.0036 

±0.071 

±0.081 

±0.031 

±4.0 

±8.8 


Table 6. Estimates of the observed standard deviations ffobs.a of the power threshold 
Pthr, expressed as pereentage errors, for the maehines eontributing to the IAEA02 dataset, 
obtained using power-law sealing. The bootstrap averages are given, as well as the 95% 
eredible intervals (Cl). 



ASDEX 

AUG 

CMOD 

DHI-D 

JET 

JFT-2M 

JT-60U 

PBXM 

Average (%) 

35.8 

21.2 

20.4 

15.9 

22.4 

15.7 

22.3 

27.7 

Cl (%) 

±9.1 

±4.3 

±3.4 

±2.4 

±3.8 

±2.2 

±4.6 

±8.1 


It should also be mentioned that, in obtaining Table 6, we again ealeulated relative errors from the 
observed standard deviations estimated by GES. However, this time, the relative errors are not the same 
for all measurements eoming from a single maehine, so we ealeulated an average for eaeh maehine 
(and similar for the eredible interval). The resulting errors on Pthr are relatively similar to those using 
log-linear sealing, with an average over all deviees of 22.7%, whieh is again higher than the 15% 
expeeted from measurement error only. 

5.3. Influence of Error Bars 

In the last eouple of experiments, we intended to assess the sensitivity of the regression analysis on 
the aeeuraey of the error bars on the ITPA data. A systematie study of this influenee is outside the seope 
of this paper, and as a first simple test, we doubled the error bars on all root ITPA variables (basieally the 
eleetron density and the magnetie field together with various geometrieal plasma parameters and sourees 
of input power), whieh were used for ealeulation of the variables involved in the power threshold sealing 
law. On average, over all maehines, this resulted in the following derived error bars: 9% on n^, 2% on 
Bt, 5% on S and 32% on Pthr- Again, these are all below the maxima quoted in Seetion 3. 
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We then performed power-law regression with MAP and GLS on the ITPA data using these larger 
error bars; the results are given in Table 7 [36]. It is observed that, based on MAP, the predietions 
for ITER are lowered relative to the analysis with the original error bars in Seetion 5.2. In eontrast, 
the predietions by GLS remain about the same as before. On the other hand, the GLS estimates of 
the observed standard deviations, listed in Table 8, are inereased for all deviees. This is how GLS 
aeeommodates the inereased error bars on the data. 

Table 7. Estimates of regression parameters and predietions for ITER in power-law sealing 
on the original seale of the H-mode threshold power using the IAEA02 dataset with all error 
bars (on the root quantities) doubled. 









Method 

Po 

|3i 

132 

133 

Pthr, 0.5 (MW) 

Pthr, 1.0 (MW) 

MAP 

0.0436 

0.581 

0.828 

0.900 

41.0 

61.3 

GLS 

0.0393 

0.725 

0.742 

0.990 

52.1 

86.2 


Table 8. Estimates of the observed standard deviations ffobsA of the power threshold Pthr, 
expressed as pereentage errors, for the maehines eontributing to the IAEA02 dataset with all 
error bars doubled, obtained using power-law sealing. 


ASDEX 

AUG 

CMOD 

DIII-D 

JET 

JFT-2M 

JT-60U 

PBXM 

49.5 

35.9 

31.7 

24.9 

32.9 

27.6 

38.9 

47.7 


In another simple test, we ehanged the error bars on fie, S and Pthr to values eomputed from 
the average pereentages mentioned earlier in Seetion 3: 4% for fie, 1% for Bt, 3% for S and 15% for 
Pthr- These are averages over all maehines, rendering the final absolute error bars (standard deviations), 
eomputed from the relative errors, less preeise. The estimation results using power-law regression with 
MAP and GLS are shown in Table 9. The results of both methods are elearly affeeted by the averaging 
step, but again, MAP is seen to be more sensitive to the ehange in the error bars eompared to GLS, whieh 
maintains estimates in a similar range as those given in Tables 3 and 5. The estimates of the observed 
standard deviations, given in Table 10, are adjusted aeeordingly by GLS. 

Table 9. Estimates of regression parameters and predietions for ITER in power-law sealing 
on the original seale of the H-mode threshold power using the IAEA02 dataset with averaged 
error bars. 









Method 

(3o 

|3i 

132 

133 

Pthr, 0.5 (MW) 

Pthr,1.0 (MW) 

MAP 

0.0488 

0.552 

0.807 

0.862 

35.1 

51.5 

GLS 

0.0429 

0.647 

0.780 

0.938 

45.7 

71.5 
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Table 10. Estimates of the observed standard deviations ffobs,® of the power threshold Pthr> 
expressed as pereentage errors, for the maehines eontributing to the IAEA02 dataset with 
averaged error bars, obtained using power-law sealing. 


ASDEX 

AUG 

CMOD 

Din-D 

JET 

JFT-2M 

JT-60U 

PBXM 

29.6 

19.1 

20.5 

19.5 

22.5 

18.1 

18.7 

20.4 


5.4. Discussion 

Several interesting observations ean be made from the experiments regarding the power threshold 
sealing in this seetion. Eirst, eonsidering Tables 3 and 5, it should be noted that there are several 
instanees where the regression parameters estimated by OES differ signifieantly from those obtained 
by GES. Eor log-linear regression, this is partieularly the ease for the dependenee of the power 
threshold on density and surfaee area and for the predieted power thresholds for ITER, as shown by 
the non-overlapping eredible intervals. Eor power-law regression, the differenee is rather situated in the 
dependenee on the magnetie field. In this ease, the power thresholds predieted by OES are also quite 
different from the results given by GES, but this time, the eredible intervals on the OES estimates are so 
wide, that they overlap with those obtained from GES. Apart from this diserepaney, the three methods 
provide eomparable absolute error bars on their estimates. 

Eurthermore, we see that the eorrespondenee between GES and MAP is signifieantly better, although 
the remaining differenees beeome particularly clear for the predicted power at higher density in ITER. 
The estimate by GES is higher than that provided by MAP, especially for power-law scaling. 

In addition, and quite remarkably, when comparing the coefficient estimates and predictions obtained 
by GES between the linear and nonlinear case, relatively consistent results can be noted. The same goes 
for the MAP estimates. In contrast, OES provides widely different results, depending on whether a linear 
(log-transformed) or nonlinear (power-law) model is used. The relatively good consistency of the GES 
estimates across regression models is a solid argument in favor of the method. 

Another noteworthy point comes from the results of the two additional tests with increased and 
averaged error bars. They indicate that for MAP (and maximum likelihood) regression, reliable estimates 
of the variability of the measurements is important. However, as discussed in Section 3, the standard error 
bars that were used in the analysis in Sections 5.1 and 5.2 are small compared to the actual variability of 
the data around the theoretical scaling law. Hence, one could speculate whether the results of the MAP 
analysis are in fact trustworthy, given its sensitivity to the error bars on the data. Therefore, at least for 
MAP, it would be better to encode the available information on the error bars in sufficiently wide prior 
distributions (which, incidentally, would be possible for GES, too). 

A related comment is that GES is clearly less vulnerable to inaccurate error specification compared to 
MAP. The mechanism behind this behavior is similar to the one that makes GES less sensitive to outliers, 
i.e., the observed standard deviation is able to capture deviations from the expected data variability with 
respect to the model. In the simple implementation of the GES method used in the present paper, the 
distinction that is made between the modeled and observed standard deviation is the main difference 
between GES and MAP. 
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6. Conclusions 

Regression and sealing laws represent erueial tools in seienee in general and in the analysis of 
eomplex physieal systems in partieular. We have presented geodesie least squares regression (GLS) 
as a method that is able to handle large uneertainties on the data and on the regression model, and 
we have demonstrated its applieation to power-law regression. Operating on a manifold of probability 
distributions, GLS has the advantage that its results ean be easily visualized in the ease of the univariate 
Gaussian distribution. However, GLS is suffieiently flexible to allow taekling mueh more general 
regression problems within the same framework. 

We have shown two examples of the enhaneed robustness of the method using synthetie data. 
GLS showed a better stability in the presenee of outliers and under a logarithmie transformation of a 
power-law, eompared to established teehniques. In addition, we have addressed the sealing of the L-H 
power threshold in magnetieally-eonfined fusion plasmas. On the basis of data from a multi-maehine 
database, it was observed that geodesie least squares provides estimates of regression parameters and 
predietions that are eonsistent aeross different regression models, in eontrast to ordinary least squares. 
Furthermore, beeause GLS allows the data uneertainty predieted by the model to be different from the 
empirieally observed uneertainty, whereas with maximum a posteriori they are, by design, the same, 
GLS is more flexible and robust at the same time. As a eonsequenee, the degrees of freedom provided 
by the parameters of the regression model better serve their aetual purpose: to parameterize a model that 
best deseribes a trend in the data, with minimal distraetion by the data “noise”. 

In future work, we intend to present a more general formulation of geodesie least squares, targeted 
at a wider elass of regression problems. In addition, various theoretieal performanee issues need to be 
addressed, ineluding uniqueness and eonvergenee properties of the optimization problem, asymptotie 
behavior of the parameter estimates, etc. On the praetieal side, we aim at establishing a broader basis for 
the performanee of the GLS method on simulated data. This should inerease the eonfidenee over a wider 
range of regression problems, as well as deviations from the regression model. 

Finally, although we have noted that GLS performs regression on a probabilistie manifold, we have 
aetually made little use of the geometrieal strueture of the manifold, save for ealeulating geodesie 
distanees. Nowadays, there are various sehemes, more sophistieated than a least-squares approaeh, 
to perform regression on manifold-valued data. From that point of view, one ean expeet advantages 
of a method performing regression between probability distributions, eaeh of them eontaining more 
information than struetureless data points in a Euelidean spaee. One possibility that we will explore in 
future work is a Bayesian regression method on a probabilistie manifold, by deseribing the distribution 
eorresponding to the regression model intrinsieally on the manifold [37]. At the same time, this will 
provide uneertainty estimates on the parameters through the posterior distribution. 
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